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Abstract 

The inelastic incompressibility is a typical feature of metal plasticity /viscoplasticity. 
Over the last decade, there has been a great amount of research related to construc- 
tion of numerical integration algorithms which exactly preserve this geometric prop- 
erty. In this paper we examine, both numerically and mathematically, the excellent 
accuracy and convergence characteristics of such geometric integrators. 

In terms of a classical model of finite viscoplasticity, we illustrate the notion 
of exponential stability of the exact solution. We show that this property enables 
the construction of effective and stable numerical algorithms, if incompressibility is 
exactly satisfied. On the other hand, if the incompressibility constraint is violated, 
spurious degrees of freedom are introduced. This results in the loss of the exponential 
stability and a dramatic deterioration of convergence behavior. 
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Nomenclature 



Ci inelastic right Cauchy-Green tensor (see fl25l) ) 

T 2nd Piola-Kirchhoff tensor (see ([27])) 

1 second-rank identity tensor 

A ■ B = AB product (composition) of two second-rank tensors 

A : B scalar product of two second-rank tensors 
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A ® B tensor product of two second-rank tensors 

II A II I2 norm of a second-rank tensor (Frobenius norm) 

II A II* induced norm of a second-rank tensor (spectral norm) (see ([1]) 

(■)° deviatoric part of a tensor 
transposition of a tensor 
inverse of transposed 

tr(-) trace of a second-rank tensor 

(■) unimodular part of a tensor (see ([2])) 

sym(-) symmetric part of a tensor 

(x) MacCauley bracket (see (fTSjl ci) 

ipe\ specific free energy density 

Dist(-, ■) "distance" between two solutions (see (139!) ) 

K yield stress 

Ai proportionality factor (inelastic multiplier) (see 

/ overstress (see (EH]) 2) 

^ norm of the driving force (see ([28]) 3) 

Sym space of symmetric second-rank tensors 

M invariant manifold (cf. ([5]), ( l30l) ) 

Pa mass density in the reference configuration 

k bulk modulus (see (133]) ) 

fi shear modulus (see (133]) ) 



1 Introduction 



The mechanical processing of materials may involve very large inelastic deformations. 
For instance, for equal channel angular extrusion of aluminum alloys, the introduced 
accumulated inelastic strain usually varies between 100 and 900 Percent (depending on 
the number of extrusions |3l])- Even larger deformations can be introduced by some 
incremental forming procedures like spin extrusion [21] (the accumulated inelastic strain 
ranges up to 1000 Percent). Due to the highly nonlinear character of the underlying 
mechanical problem, a correct numerical simulation of such "long" processes is by no 
means a trivial task. It is desirable to have numerical algorithms which would be stable 
with respect to numerical errors, even if working with big time intervals and big time 
steps. 

The assumption of exact inelastic incompressibility is widely implemented for construction 
of material models of metal plasticity and creep (see, for instance, [H]). Extensive studies 
were carried out concerning the construction of numerical integration algorithms which 
exactly preserve the incompressibility of the inelastic flow [Hll0lll3f20ll23ll27f30ll31j P1 

^ The incompressibility condition is given by a linear invariant in the case of infinitesimal strains 
inelasticity. Since the hnear invariants are exactly conserved by most of integration procedures 
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In this paper, we asses those factors that result in a more accurate computations, especially 
when integrating with big time steps and for long times. To this end, we analyze the 
structural properties of the inelastic flow governed by a classical material model of finite 
viscoplasticity. The material model is based on the multiplicative decomposition of the 
deformation gradient into inelastic and elastic parts. For simplicity, no hardening behavior 
is considered in this paper. However, the proposed methodology can be generalized to cover 
more complicated material behavior as wellH] 

We pay especial attention to the exponential stability of the inelastic flow, which is the key 
notion of the current study. We say that the solution to a Cauchy problem is exponentially 
stable, if for small perturbations of initial data, an exponential decay estimate holds (see 
Section 2.1). From mechanical standpoint, the exponential stability implies fading memory 
behavior! ^ I Moreover, the exponential stability is deeply connected to contractivity (B- 
stability) of the system of equations, which can be used for stability analysis of numerical 
algorithms (see the monograph by Simo and Hughes [29]). 




The main conclusions of this paper regarding the problem of finite viscoplasticity are as 
follows. 

• The exact solution is exponentially stable with respect to small perturbations of initial 
data, if the incompressibility constraint is not violated. 

• In the case of exponential stability, the numerical error is uniformly bounded. In partic- 
ular, there is no error accumulation even within large time periods. 

• If the incompressibility constraint is violated by some numerical algorithm, then, in 
general, the numerical error tends to accumulate over time. 

There exists a rich mathematical literature dealing with existence, uniqueness, regularity, 
and asymptotic behavior of solutions for certain plascticity/viscoplasticity problems in the 
context of infinitesimal strains (see p!|5ffTi] and references therein). A class of material 
models of monotone type which includes the class of generalized standard materials was 
defined and analyzed in fT]. In the context of finite viscoplasticity, however, only few 
theoretical works exist. Some preliminary investigations have been made by Neff in [24j . 




In this paper, we analyze the well-known material model of finite viscoplasticity. The 
stability is proved analogously to the classical Lyapunov approach, based on the use of 
Lyapunov-candidate-functions. In fact, the hyperelastic potential is used to construct a 



(cf. [7]), the problem of the conservation of incompressibility only appears when working with 
finite strains. 

^ Using a series of numerical tests, it was shown in [27J that the use of geometric integrators 
allows to eliminate the error accumulation even in the case of a more complex material be- 
haviour with nonlinear isotropic and kinematic hardening. In general, however, the construction 
of consistent integration procedures for the finite strain inelasticity is still an open problem (cf. 



^ As Truesdell and Noll [32j put it, "Deformations that occurred in the distant past should have 
less influence in determining the present stress than those that occurred in the recent past" . 





[33]). 
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suitable Lyapunov candidate (cf. [29]). 



The paper is organized as follows. In Section 2, we define the notion of exponential sta- 
bility and prove the main theorem, which states that the numerical error is uniformly 
bounded if the exact solution is exponentially stable. A simple one-dimensional exam- 
ple is presented. In the next section, a classical material model of finite viscoplasticity 
is formulated in the reference configuration. The change of the reference configuration is 
likewise discussed. Section 4 contains the definition and analysis of the distance between 
two solutions in terms of energy (Lyapunov candidate). Next, the time-evolution of the 
distance is evaluated and the exponential stability of the exact solution is proved. Finally, 
the results of numerical tests are presented, which illustrate the excellent accuracy and 
convergence characteristics of geometric integrators. 

We conclude this introduction with a few words regarding notation. Expression a := b 
means a is defined to be another name for b. Throughout this article, bold-faced symbols 
denote first- and second-rank tensors in M'^. A coordinate- free tensor setting is used in 
this paper (cf. |15|[28]). The scalar product of two second rank tensors is defined by 
A : B = tr(AB'^). This scalar product gives rise to the norm by ||A|| := v^A^A. 
Moreover, we denote by || ■ ||* the induced norm of a tensor 



||A||* := max ||Ax||2, ||x||2 := a/x ■ x. (1) 

l|x||2 = l 

The overline (■) stands for the unimodular part of a tensor 

A := (detA)-^/='A. (2) 

The deviatoric part of a tensor is defined as A° := A — |tr(A)l. The notation O 
stands for "Big-0" Landau symbol: f{x) = 0{g{x)) clS X — ^ Xq iff there exists C < 
oo such that ||/(x)|| < C||(?(x)|| as x — ^ xq- The inequality f{x) < 0{g{x)) is understood 
as follows: there exists f{x) = 0{g{x)) such that f{x) < f{x). 



2 Differential equations on manifolds and exponential stability 

2.1 General definitions 

Let us consider the Cauchy problem for a smooth function y{t) G M" 

y{t) = f{yit),d{t)), y{to)=yo. (3) 
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Here, the initial value yo and the function d{t) are supposed to be givenQ Denote the 
exact solution to ([3]) by y{t,yQ,to). In particular, we have 

yito,yo,to) = yo- (4) 

Suppose that all solutions lie on some manifold M C 

y{t, 1/0, to) e M, for all t > to, y^ e M. (5) 

Then we say that ([3]) is a differential equation on the manifold M (cf. [6|7] ). 

Next, we say that the solution y{t) to the problem ([3]) is locally exponentially stable on 
M, if there exist 6 > 0, 7 > 0, Ci < 00, such that the following decay estimate holds 

\\y{t,yl,'\to)-yit,yi'\to)\\ < C, e'^l^ - ^0) \\y^^^ - yf\l (6) 

for all to > 0, y'i\y''o^ G M such that \\y'i^ - y{to)\\ < 5, - y{to)\\ < 6. 

We note that somewhat different interpretation of the exponential stability can be met in 
the literature as well (cf., for example. Section 2.5 of [16]). 

Next, let us consider a numerical algorithm which solves (I3l) on the time interval [0,T]. 
Denote by ""y the numerical solutions at time instances "t, where = ^t < ^t < '^t < ... < 
^t = T, and ^y = yo. Suppose that the error on the step is bounded by the second power 
of the step size. More precisely 

"y, "t) - "+^1/11 < Csr+H - (7) 

where C2 < 00 (cf. figure [1]). For simplicity, we will consider constant time-steps only: 
At = - "t = const. 

2.2 Main theorem 

With definitions from previous section we formulate the following theorem. 
Theorem 1. 

Let y{t) = y(t,yo,0) be the exact solution. Suppose that conditions ([6]) and ([7]) hold. 
Moreover, suppose that the numerical solution of problem ([3]) lies exactly on M. Then 
there exist a constant C < 00 such that 

ry - yCt) II < CAt, as At ^ 0. (8) 

Here, the constant C does not depend on the size of the time interval [0,T]. 

^ The system ([3]) is a system with input, and d{t) is interpreted as a forcing function. 
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X Numerical solution 
• Exact solution 




yi't) 



H H H H H 

Fig. 1. Analysis of error accumulation. 

Proof. The proof is a modification of the standard error analysis (cf. [2]). In this paper 
we prove the theorem under assumption that 5 = oo. The proof can be easily generalized 
to cover arbitrary values of 5 > by using mathematical induction and by assuming 
At < 7(5/(2 max((:7i,l) C2). 

First, note that = Thus, (cf. figure [I]) 

ry~y{^t)\\ < ry-yi%--'y,--H)\\ + ||y(% ""'t) - y(% ""V ""'t) II + - 

+ ||y("t,V't)-y("t,V°t)||. (9) 

Next, from ([6]) we obtain for all A; = 1, 2, n — 1 

\\yCt,''y,H)-yi%'-'y,''^h)\\<C, e-7("^ " \\yi%'y,H) -y{H,'-'y,''-h)\\. (10) 
Substituting ([lU]) in ([2]), we get 

ry-yCt)\\< 

ry-y{X"-'y,''-'t)\\ + C,j:e-^i ^) \\y{%'^y,H) - y{\'-'y,'-H)\\. (11) 

fc=i 

Obviously, y{^t, ^y, H) = ^y. Without loss of generality, we can assume that Ci > 1. Next, 
substituting error estimation ([7j) into ffTTl) . we get 



y-y 



-t)||<CiC2(At)2^e-T'(^- t). 



(12) 



k=l 



But, "'t — = {n — k)At. Thus, taking into account the well-known expression for an 
infinite geometric series {J2iZo^^ — ^/i^ ~ ^) f*^^ kl < 1 ), we get for small At 



^e-7r^-'^)<f;e-^7At^ 1 

k=l i=0 1 — e 



< 2- 



7At + 0((At)') - 7At' 



(13) 
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Fig. 2. Rheological model (a), and inelastic flow under monotonic loading (b). 
Finally, it follows from ^ 

in c 

iri/-l/("t)|| < ' ' At, asAt^O ■ (14) 



Remark 1. The proof is essentially based on the assumption that ''y G M. In general, if 
the numerical solution '^y leaves the manifold M, the decay estimation ([6]) is not valid. 

Remark 2. The theorem states that the error is uniformly bounded in the case of expo- 
nential stability. Thus, there is no error accumulation in the sense that the constant C in 
(IH]) does not depend on the overall time T. Moreover, let e > be some small value. By 
choosing At < 7e/(2CiC2) the numerical error ||"?/ — y("t)|| is guaranteed to be less than 
e. 

Remark 3. If the exponential stability is replaced by the assumption that the right-hand 
side of ([3]) is a smooth function of y, a weaker error estimation is valid (cf. [2]) 

ry - yC't^ < C e^^ T At, as At ^ 0, (15) 

where L = sup The effect of growing multiplier on the right hand side of (ITSl) is 
referenced to as an effect of error accumulation. In that case, in order to guaranty a 
sufficient accuracy, the upper bound for At must depend on T. That makes the practical 
solution of some problems extremely expensive for large values of T. 



2.3 One- dimensional example 



Let us consider a simple example which illustrates the notion of exponential stability. We 
examine the response of a one-dimensional viscoplastic device shown in Figure [2] (a). The 
closed system of (constitutive) equations is as follows: 

The total strain is decomposed into elastic part Eq, and inelastic part Ei 

E = Ee + E[. (16) 
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The stress o on the elastic spring is governed by elasticity law (i? > 0). 

a = Ee,. (17) 
The time derivative of the inelastic strain is given by 

^i = ^(/)|^, f:=\a\-K, (x) :=max(x,0), (18) 

where material constants K > Q and > are referred to as yield stress and viscosity, 
respectively. 

In order to use the results of previous subsections, we rewrite the problem in the form 

e, = e,{eue{t)) = -{E\e{t) - e^l ~ K) sign(£(t) - £;). (19) 
V 

Let el^\t) and e['^\t) be to two solutions to (fT9|) . Following [29j, we recall that 



defines an energy norm which is the natural norm for the problem under 
consideration[£j Next, we consider a monotonic loading 

e{t) =et, e = const > 0. (20) 

Let us show that the exact solution satisfying the initial condition = is exponentially 
stablelZl Without loss of generality, we can assume that to = in estimation If 
ki'^'*(0) ~ 0| ^ ^ k G {1,2}, then there exists time instance t' = t'{6) such that the 
condition / > /o > holds for both solutions {e[^\t) and el^\t)), if t > t' (see Figure [2] 
(b)). Then, under that assumption 

^%^ = -^. .,(4"..W)-^,(.P'..W) = --(4"-.P'). (21) 

osi rj rj 

Therefore, we get from fl2T|) 

(1e(.;^) - eP)^ ■ = Eiel'^ - .P))(4^) - #^) = -^(.« - ef^ f, for t > t'. (22) 



Due to the contractivity (for details see p9]), lE{e[^\t')-e'f\t')y < ^E{e[^\o)-e['^\o)Y. 
Moreover, integrating fl22|) over [t',t] and taking the contractivity into account, we get 

lEiel'\t)-ef\t)r < ii^(.P^(0) -ep)(0))^e-f (23) 



5 It is known (see [29]) that ^E{e\ '{t) - e[ '{t) f is not increasine. This effect is referenced to 
as contractivity. 

^ For the current example, the geometric property ?/ G M is trivial: we put M = M. 
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Taking the square root of both sides we obtain the required exponential decay estimation 
Et' p 
with Ci = < cx) and 7 = ^ > 0. 



3 Material model of multiplicative viscoplasticity 

Let us consider a classical material model of finite viscoplasticity (see, for example, 
3.1 Constitutive equations 

The model is based on the multiplicative split of the deformation gradient F 

F = FcFi. (24) 

Here, Fe and Fi stand for elastic and inelastic parts, respectively (see [TTfTS] ). The mul- 
tiplicative split can be motivated by the idea of a local elastic unloading. A somewhat 
more consistent motivation can be derived from the concept of material isomorphism [3]. 

Along with the well-known right Cauchy-Green tensor C = F^F, we introduce a strain- 
like internal variable (inelastic right Cauchy-Green tensor) as 

Ci = F^Fi. (25) 



In this paper we consider strain-driven processes. More precisely, we assume the defor- 
mation history C(t) to be given. The material response in the time interval t G [0,T] is 
governed by the following ordinary differential equation with initial condition 



Ci = 2|(CT)°Q, Ci|t=o = C°, detC? = 1, C° G Sym. 



(26) 



Here, the 2nd Piola-Kirchhoff tensor T, the norm of the driving force ^, and the inelastic 
multiplier Ai are functions of (C, Ci), given by 



2p, 



dC 



Ci=const ' 



(27) 



f:=d 



tr 



~ \ D 



(CT) 



(2S 



The material parameters > 0, ?7 > 0, m > 1, K > 0, and the isotropic real-valued 
function ip^^i are assumed to be known; /cq > is used to get a dimensionless term in the 
bracket. 
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Remark. The right Cauchy strain tensor C is symmetric. Since the function ifjei is 
isotropic, it makes no difference whether the derivative in (127|) is understood as a general 
derivative or as a derivatives with respect to a symmetric tensor (cf. [28j ) . 

Next, we remark that the right-hand side in (!26|) ^ is symmetric (cf. [28]). Moreover, 
taking into account the property tr(AB) = tr(BA) and combining the Jacobi formulae 

^'^^A^'^ ~ det(A) A^""^ with the evolution equation (12^ -^. we get 

(detQ)- = : Q = 2|det(Q) Cr^ : (ct)°Q = 2|tr[(ct)°] = 0. (29) 



Therefore, the exact solution of (]26l) - (12811 has the following geometric property 

Ci G M, M := {b G Sym \ detB = l}. (30) 

We note that the current material model is thermodynamically consistent. That means 
that the Clausius-Duhem inequality holds for arbitrary mechanical loadings 

6, := -^t : C - (^ei(CCi-i))' > 0. (31) 

In particular, we get a reduced inequality for relaxation processes (C = const) 

[CCr'))' <0. (32) 
One mathematical interpretation of this inequality will be discussed in Section 4.1. 

To be definite, we use the following expression for the free energy density ipei (generalized 
Neo-Hooke model [IT] ) 

PnM^) ■= ^{\nVd^y + |(trA - 3), (33) 

where k > 0, fi > are known material constants (bulk modulus and shear modulus, 
respectively). 

Substituting ([3SD in (EZD we get the 2nd Piola-Kirchhoff stress tensor in the form 



T = A; In^det(C) C^^ + /i C-\CCrY- (34) 
In what follows we analyze the exponential stability of the exact solution Ci{t). 
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3.2 Change of reference configuration 



In order to simplify the analysis of the material model, we may need to rewrite the 
constitutive equation with respect to some "new" local reference configuration Fq. In 
what follows, we suppose that this configuration is isochoric, i.e. det(Fo) = 1. The "new" 
deformation gradient, right Cauchy tensor, and inelastic right Caushy tensor are given by 

pnew.^pp-i^ C''^" := Fq ^CFq \ Cf"" := Fq TQFo \ (35) 

The 2nd Piola-Kirchhoff tensor T, the norm of the driving force the inelastic multiplier 
Ai, and the overstress / are transformed as follows 

,j,new _ p^^pT^ jncw_^^ A^" := A; (36) 

Since ifjei is isotropic, ipexi^AS) = '?/'ei(BA). Using that property, it can be checked that 
ipei{CCi~^) is invariant under the change of reference configuration 

^^^^C--(Cr")-') = ^ei(CCr'). (37) 

The closed system of equations with respect to the new reference configuration is obtained 
from (|26ll — ( |28ll by replacing all quantities by their "new" counterparts. 



4 Analysis of exponential stability for multiplicative viscoplasticity 

4-1 Measuring the distance between solutions in terms of energy 

Suppose that Cp^ (t) and Cp^ (t) are two solutions to the problem fl26l) — fl28l) (with the 
same forcing function C(t)). Next, suppose that there exists a constant L < oo such that 

\\{Ciy/\t)\\ < L, ||(Cpyi/2(t)|| < L for alH > 0, e {1,2}. (38) 
We introduce the following measure of distance between two solutions in terms of energy 

Dist(C;^), Cp)) := yp,^,i(C«(CpVi). (39) 
This measure has the following properties: 



' The relation (|39|) can be seen as a generalization of the energy norm A/iS(ep^ -ep^)2 (cf. 
Section 2.3). 
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(i) Invariance under the change of reference configuration 



Dist (C 



new ^^(2)Njncw 



Dist(Cp\Cp^ 



(40) 



ii) For small Cp'' — Cp"*, there exist constants Cci > and C4 < oo such that 



C3 ||Cp^ - Cp^ll < Dist(q'^ C\'>) < C4C\'> - C\ 



(1) n(2)^ 



^(1) n(2)i 



(41) 



(iii) For all C[^\t), Cp^(t) e M we have Dist(cS^\ Cp^) > and 
Dist(cS^\ CS^^) = 0, if and only if C[^^ = 



(42) 



Proof. 

(i) : Identity fj40|) can be proved similarly to the invariance property fl371) . 

(ii) : First, it follows from fl33p that for small A we have (see Appendix A) 

Pa^ei(l + A) = ^(trA)2 + f tr((A°)2) + ©(A^), 



(43) 



where A = ||A||. Note that tr((A°)2) = ||A°f for A G Sym. Thus, |(trA)2 + 

^tr(^(A^)^^ is a norm on Sym. Since all norms on Sym are equivalent, there exist con- 
stants C3 > 0, C4 < cxD such that for small A G Sym we have 



C^||A||2<p^^,i(l + A)<C;||Af. 



Next, due to the property 
we have 



^el(AB) = 7/>el(BA) 



Dist(cP\ cp)) = ^pH^ei((CpV^/' (CpV^/^ 

Moreover, taking into account that ||AB||^, < ||A||^, ||B||*, and that the norms 
II ■ II are equivalent, we get 

||ABC|| < CIIAII ||B|| ||C||, 
with some constant C < 00. Thus, 



(44) 

(45) 

(46) 
and 

(47) 



cp))-i/2 cp) fcp))-i/2 _ 1 ^ rcp))-v2 fcp) - cFl (cPh-^/^ 



1 V ~ 1 



103 , 
< C 



.Cp))-l/2 



cp) - cp) 



(48) 
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< C ' 



CpV^' ' (CS'V'/' (CpV^/^ - 1 . (49) 



Further, substituting A = (C[^V^^^ C[^^ (C[^V^^^ - 1 in (jS]), and combining it with 
), and (HHD we get 



C^\\iC[y/'\\-' ||Cp) - Cp^ll < Dist(cP, CP^) < C4{C['Y'^'\\' ||CP - Cp)||. (50) 
Finally, combining ( l50l) with ( l38l) we obtain fHTi) . 

(iii): We note that '?/'ci(A) > 0. Moreover, i/je\{A) = if and only if A = 1 ■ 

In view of properties (i) — (in), the function Dist is a natural measure of distance for the 
problem under consideration] ^ I 

Moreover, the dissipation inequality (l32i) . which holds for all relaxation processes, can 
be interpreted as follows: during relaxation, the distance (measured in terms of energy) 
between any solution Cp"* and a constant solution C[^^ = 1 is not increasing. 

4-2 Sufficient condition for exponential stability 



Let us consider a loading program (strain-driven process) {C}teio,T] on the time interval 
[0,r]. Let Cp\ C[^^ e M be two solutions. In order to prove the exponential stability, it 
is sufficient to prove that there exists t' > and 7 > such that for all t > t' (cf. (1221) ) 



Dist(Cp\Cp))2)' < -7 Dist(CP\Cl 



(1) r^{2)N2 



(51) 



Indeed, in that case, using the Gronwall's inequality we get from fl5T|) the following decay 
estimation 

Dist(C['^(t), Cp^(t))2 < Dist(C['^(t'), Cp^(t'))' e-^(*-*') . (52) 
Combining this result with fHTj) . we get the required estimation of type ([6]). Thus, the 
uniform error estimation of Theorem 1 follows immediately from fl5^ . 



4.3 Reduction of the stability analysis to a simplified problem with C = 1 



Let t'^ be an arbitrary time instance. In this section we discuss a procedure, which helps 
to simplify the examination of the inequality (ISTl) at time t^. 



^ The function Dist is not symmetric: Dist(A,B) 7^ Dist(B,A). Symmetrized functions 
ca n be defined by Di stf^'"(A, B) := l/2(Dist(A, B) + Dist(B,A)), Dist2^"'(A, B) := 
Y/Dist(A, B)Dist(B, A). Nevertheless, none of these functions determine a metric on M, since 
the triangle inequality does not hold. 
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The first simplification of the problem is as follows. We note that quantities 
Dist(Cp^(tO),Cp^(tO))2, and (Dist(Cp^(tO), Cp)(t0))2)' ^^p^^^^ g^j^iy ^^o)^ Cp^(t°), 

and Cp'*(t°) but not on C{t^)- Therefore, at the examination of fl5T]) at t = t*^ we can 
replace the actual loading programm {C}fg[o,T] by a constant loading (relaxation process): 
we take a constant C(t°) instead of loading C{t), where (■) stands for a unimodular part 
of a tensor. 

The second simplification is as follows. Let Fq be some "new" reference configuration and 
det(Fo) = 1. There is a one to one correspondence between the solutions C(^\t), Cp'*(t) 
of the problem with the forcing function C{t) to the solutions (CS^^)''<=™(t), (C[^^)"°*(t) 
with the forcing function ^^^"(t) (cf. Section 3.2) 

C"^™(t) = Fo-TCFo 1, (Cp))"^-(t) =FoTCp^(t)Fo\ ke{l,2}. (53) 

It follows from that 

Dist((Cp^)'^^"(t°), (Cp^)"'="(t°)) = Dist(Cp^(t°), Cf\f)), (54) 



Dist((Cp))-"(t), {Cl'Y-{t)y]\,=,o = [Dist(cp)(t), Cp)(t))'] |,=,o, (55) 
Therefore, estimation flSTl) is equivalent to 

Dist((Cp^)'^^"(t), (Cp^)'^^"(t))'] |i=iO < -7 Dist((Cp^)'^^"(t°), (Cp^)°^"(t°))'. (56) 

1 /2 

Without loss of generality we assume det(C(t°)) = 1. By choosing Fq = (C(t°)) the 
problem can be reduced to the simplified problem with C{t^) = 10 



44 Evaluation of (m(cS^\ CfV)' 



In this section we evaluate (Dist(C[^\Cp^)2)' at some fixed time instance t°. Without 
loss of generality (cf. the previous section) it can be assumed that C(t°) = 1. In that 
reduced case, the evolution equation (126!) takes the form 



Q = a(||(Crni)(Cr^)°Q, aix):=^( ^'' f'^^ Y ■ (57) 



^ Alternatively, the problem can be reduced to the case c[^\t^) = 1 by choosing Fq 
(l)^.0^^1/2 
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Next, using the product rule we get from ([5 



(cp)cp)-i 



^{l)^(2)-l , ^{l).^(2)-l^ 



(^{l)Qa)-lQ{l)Q{2)-l _^ (^(l)^(^(2)-ly^(2)^(2)-l 



«(i)(Cp)-^)°-a(2)(cp)-i)D 



i(2)-l 



(5^ 



where a(^) := a(||(Cp^"')°||) for A; G {1,2} 



Further, we compute the derivative of V^ei(A) using a coordinate- free tensor setting (see, 
for example, [ISII2B])- 



P»^ = |lnV«AA-+^A-(AY. (59) 

We abbreviate A := ||CS^^ - Cp^||. Note that (see Appendix B), since C['^\ G M 



tr((Cp)-^-Cp)-^)C«) = 0(A2), tr(c«(Cp)-^-Cp)-^)) = 0(A2),as A ^ 0. (60) 
Thus, using ( l60!) i we get 



Cp)-ic«l° = ((Cp)-^ - Cp^-^)C« + l)° = (Cp^-^ - Cp)-^)Cp) + 0(A2). (61) 



1 ^1 



Combining ( l59l) with ( ]6T|) we get 



a^ei(c;^)cp)-^) 
^ 9(cp)cp)-^) 



(^(l)(-;.(2)-l 



D 



= ^Cp)-^Cp)(Cp)-^ - C«-^)C« + 0(A2) = |(Cp)-^ - Cp)-^)CP) + 0(A2). (62) 
Next, denote by d the derivative of a{x) at X = ||(Cp^"^)°||. Therefore, 



a« - a(^) = d (||(C«-^)°|| - 11(0^-^)^11) + O(A^) 



(i)-i^D . fn{i)-i _ c(2)-i 



II (Q 



(1)- 



0(A2). (63) 



It can be assumed that the overstress / = yu||(CP'' '^)^|| — is bounded by \I\K. 

Thus, we suppose ^J^Kj < IKCp'' "'^)^|| < 2^K/fi. Here, the first inequality is needed 
to ensure the overstress is larger than zero. 
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Using the property A : (BCD) = (B'^AD'^) : C it follows from ([58]) and ([62]) that 

rj , /r^{l)r^(2)-l\ 

(D,st(cp'.cp')7 » (m4.(c!"cP'-'))- = P. ^!;:g; J;;,. ' : (cfcp--)- 

•^(^i ) 



2 



J^(i)(Cp)(cp)-i - Cp)-^)) : (CP)-^ - Cp 



_ - a(2))(cp)(cp)-^ - CP)-^)) : Cp)-^ + 0{A') ^Fi + Fjj + 0{A'), (64) 

where Fj and F// are given by 

F, := -|„mtr((cr'-' - CP'-')C<"(C™-' - CP'-')), (65) 

^" -f i.wthoii K''')" ^ (C,"'- - Cp'-))(l : (Cp'-' - cP'-')). (66) 

IK^i J II 

Now, for any pair of real positive numbers {6, A) let us define a subset of M x M by 

S{e,A) := {(Cp\Cp^) G M X M I ||(Cp^"')°|| < e, ||Cp^ - Cp^|| < A}. (67) 



By definition, put 



dF, tr((cp)-^ - Cp)-^)Cp)(C«-^ - Cp)-^; 

(68) 

There exists a function 5(6*) > such that 

q{e) > <l>(Cp\ Cp^) + 0(A) for all (Cp\ Cp^) G 5(^, A). (69) 



The numerical evaluation of the function q{9) is discussed in the Appendix C. Moreover, 
suppose that 

> q(2J-K/ij)d. (70) 



. 3 

This condition will be discussed in the next section. Multiplying both sides of flTUl) 
by Fi^ < and noting that ^FjO(A) = 0{A^) we get for all (Cp\Cp^) G 
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S(2j|A7;i,A) 



F, f ,(2f-K/,)j^F, ? (*(C|". CP') + 0(A)) A^F, S -2F„ + 0(A= 

Multiplying both sides of flTTl) by 1/2 and adding 1/2F/ + F/j, we get 

+ < 1/2F, + 0(A3). 
Combining this result with (IMI) we obtain 

(Dist(Cp\ Cp^)2)' < 1/2F, + 0(A3). 
Next, if / > /o for some /o > 0, then there exists C5 > such that 



2 



(CI 



(i)-i 



Cp)-i)(cp))i/2 



Therefore, for small A, inequality f[73l) yields 

(Dist(Cp\cS'^)2)' < 1/4F,. 
Similarly to the proof of ( HTl) we obtain with some Ce > 

(ca)-i_cp)-i)(cp))i/2 ' 



< 



= -f ««(Ca/C:)||(C«)Vi'Wi(cPcp)-p 



(71) 



(72) 



(73) 



(74) 



(75) 



(76) 



Finally, combining (1751) with ( 1761) we get the required estimation (|5T1) if the following 
assumptions hold: < /o < / < \/f-f^, a^^^ > g( 2A/|i^'//i Id. 



^.5 Analysis of the sufficient stability condition 



In this section we analyze the condition (ITOll which was used in the previous section to 
prove the inequality (ISTI) . First, we suppose ||(Cj~^)°|| > ^J^Kj ^ to ensure the overstress 
is larger than zero. Using fl57|) 9 it can be easily shown that fITOl) is equivalent to 

||(Cr^)°|| >a;e., (77) 
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No flow Exponential stability / 
•— • • ► 



3 



0/c, 



3^^ 



Fig. 3. Domain of exponential stability, 
where the critical value Xcr is given by 



Ik + (m - l)fj.q{9) + ^ {^K + (m - 1)^(9)^ + 4/xg(^)^ir) /(2/i) , (78) 



with 6 = 2J^K/ ji. For small values of q{2Ji^K/ fx)^ a simple estimation for Xcr is valid 



Xr 



'-K/f^ + mg(2^ir//i) + 0((g(2^ir//i))2) . 



Alternatively, in terms of the overstress /, the condition (1701) is equivalent to 

/ ^ fcri 

where the critical overstress fcr is estimated by 



(79) 



^0) 



= m/ig(2y|ir//i) + 0((g(2y|ir//i))- 



^1) 



The situation is summarized in figure [31 

For instance, for aluminium alloy we put K = 300 MPa, fi = 25000 MPa. Thus 2^J^K/ fj, ^ 
0.014. Next, g(0.014) ^ 0.00000023 (See Appendix C). Therefore, the critical overstress 
is given by fcr ~ 0.0057 MPa. For physically reasonable values of m (m < 100) this 
critical value is negligible compared to the size of the elastic domain \f^K ^ 245 MPa. 

Remark. Since the overstress / is isolated from zero due to the sufficient stability con- 
dition (IHOl) . the current theory can not be applied to exactly quasistatic processes. On 
the other hand, the theory is directly applicable for nearly quasistatic processes with the 
oversress / larger that fcr- 



5 Accuracy testing of implicit integrators 



The numerical implementation of the material model (p6|) — (128|) within a displacement 
based Finite Element Method (FEM) with implicit time stepping is based on the implicit 
integration of the evolution equation (1261) (see, for example, [29]). This procedure should 
provide the stresses as a function of the strain history. 
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More precisely, suppose that the right Cauchy-Green tensor "+^0 at the time tn+i = 
tn + At is known and assume that the internal variable Ci at the time t„ is given by "Ci. 
We need to compute the internal variable Ci at the time in order to evaluate the 
stress tensor "+it = T("+iC, "+^0;). 

Note that the norm of the driving force and the overstress / can be represented as 
functions of "+^C and "+^0;: 



5("+iC,"+^Ci) = ytr[("+iC t("+iC,"+iCi 



/("+^C,"+^Ci) = 5("+^C,"+^Ci) - \ -K. 



12) 
i3) 



For what follows it is useful to introduce the incremental inelastic parameter 

e := At "+^Ai. 



^4) 



Thus, according to the Perzyna rule, we get the following equation with respect to C, 
"+iCi and C 

\ iT^) ■ ^85) 

The remaining equation for finding unknown ""'"^Ci and ^ is obtained through the time 
discretization of (!26l) . which will be discussed in the next section. 



, At/ /r+^c, 



5. 1 Euler Backward Method and geometric implicit integrators 



We introduce a nonlinear operator B("+^C, "+^Ci, ^) as 



(86) 



Let us consider the classical Euler-Backward method (EBM) (see, for example, [IfQpQ] ) 
being applied to the evolution problem (!26l) 



"+^Ci = [l - B("+iC, "+^Ci, O] "Ci. (87) 

Since the symmetry of the internal variable ""*"^Ci is exactly preserved by the EBlvf^. 
this equation is equivalent to 

"+iCi = sym( [l - B("+^C, ^^+^0;, O] "Q) . (88) 



Moreover, it was shown in |27j that the symmetry is exactly preserved by Euler-Backward 
method and Exponential Method even in a more general case of a nonlinear kinematic hardening. 



19 



The modified Euler-Backward method (MEBM) (see [T3|27] ) uses the following equation 



"+iCi = sym([l - B("+iC,"+iCi,0] ^ "Q 



^9) 



Finally, the Exponential Method (EM) (see, for instance, [4|22f23f35j ) is based on the 
use of the tensor exponential exp(-). As it was shown in [27], the Exponential Method can 
be written in the following form: 



"+iCi = sym(exp [b("+iC, "+iCi, O] "Q). 



(90) 



Combining fl85|) with one of the discretization methods (equations flHHl) . fl89|) or fl90|) ) a 
closed system of equations is obtained. One possible solution strategy for the resulting 
problem was discussed in [27j, and the application of a coordinate-free tensor formalism 
to the numerical solution was analyzed in [28]. 

We note that the geometric property of the exact flow (Ci G M) is exactly satisfied by 
MEBM and EM. Therefore we refer to these two methods as to geometric integrators. On 
the other hand, the incompressibility constraint is violated by the classical EBM. 

For all the three methods, the error on the step is bounded by the second power of the 
step size (cf. estimation ([7])), if the right-hand side is a smooth function. Strong local 
nonlinearities due to the distinction into elastic and inelastic material behavior or due to 
the non-smoothness of the loading function C{t) may increase the error on the step. 



5.2 Testing results 



The theoretical results obtained in this study are validated via a series of numerical tests. 
Let us simulate the material behavior under strain controlled, nonproportional and non- 
monotonic loading in the time interval t G [0, 300]. Suppose that the deformation gradient 
is defined by 

F{t) = F(t), (91) 

where F'{t) is a piecewise linear function of time t such that F'(0) = Fi, F'(IOO) = F2, 
F'(200) = F3, and F'(300) = F4 with 



Fi := 1, F2 : = 










1 


V2 








\ 



1 

V2J 



^lo^ 



1 
1 



2 
1 

0^, 
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Fig. 4. Accuracy analysis concerning Euler Backward Method (EBM), Modified Euler Backward 
Method (MEBM), and Exponential Method (EM). 

Thus, we put 

r (1 - t/100)Fi + (t/100)F2 if t e [0, 100] 
F'(t) := J(2-t/100)F2 + (t/100-l)F3 if t G (100, 200] . 
[ (3 - t/100)F3 + (t/100 - 2)F4 if t G (200, 300] 

The material parameters used in simulations are summarized in table [H 
Table 1 

Material parameters 
k [MPa] 1^1 [MPa] K [MPa] m [-] r] [s^^] ko [Mpa] 

73500 28200 270 3.6 2-10^ 1 

Next, we suppose that the reference configuration is stress free. Therefore we put 

Ci|i=o = 1. (92) 

The numerical solution obtained with extremely small time step (At = 0.01s) will be 
named the exact solution and denoted by Cf^"'^*. Next, the numerical solutions with 
At = Is and At = 0.5s are denoted by C™''. The error ||C™'' - Cf^'""*|| is plotted on 
figure m 

For all three methods the error is proportional to At. Moreover, in accordance with 
Theorem 1 (cf. Section 2.2), the error is uniformly bounded for geometric integrators 
(MEBM and EM). More precisely, the error is bounded by CAt, where the constant C 
does not depend on the size of the entire time interval. Next, since the incompressibility 
condition is violated by EBM, the geometric property (!30!) is lost and some spurious 
degrees of freedom are introduced. In that case, only a weaker error estimation is valid: 
II (~<nMmer_(~<exact 1 1 ^ C'(T)At, whcrc C{T) dcpeuds ou the size T of the entire time interval. 
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incompressibility 
is satisfied 



solution is 
exponentially stable 



the numerical error 
is uniformly bounded 



no error accumulation 



incompressibility 
is violated 



spurious degrees of 
freedom are introduced 



exponential 
stability is lost 



the numerical error tends 
to accumulate over time 



Fig. 5. Summary of the main results. 
6 Discussion and conclusion 



In the last decade, intensive research has been carried out concerning the development of 
so-called geometric integrators for the evolution equations of finite plasticity/ viscoplasticity, 
which exactly preserve the inelastic incompressibility condition. The excellent accuracy 
and convergence properties of such algorithms were analyzed by numerical computations. 
Particularly, the long term accuracy of geometric integrators was analyzed in the paper 
[27] , and the absence of error accumulation was numerically verified. In the current study, 
a rigorous mathematical formulation of this phenomena is proposed. The main result of 
the current paper is as follows: the numerical error is uniformly bounded by CAt if the in- 
compressibility condition is satisfied. In terms of a classical model of finite viscoplasticity 
we prove that all first order accurate geometric integrators are equivalent in that sense. 
This theoretical result corresponds with the numerical tests. Indeed, MEBM and EM are 
equivalent concerning the accuracy and convergence (cf. figure H]). The main results are 
summarized diagrammatically on figure O 

The property of the exponential stability of the exact plastic flow was mathematically 
analyzed in this paper. Obviously, that property must be utilized during the development 
of new material models and corresponding algorithms in order to improve the accuracy 
and convergence of numerical computations. 



Appendix A 

Suppose A = II A II ^0. Let us show that 

PH,^ei(l + A) = ^(trA)^ + f tr((AD)2) + ©(A^). (93) 

First, recall the Taylor expansion of det(l + A) up to second order 

det(l + A) = 1 + tr(A) + l/2(tr(A))2 - l/2tr(A2) + 0{A^). (94) 

Therefore, 

^det(l + A) = 1 + l/2tr(A) + ©(A^), (95) 
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-(lnv'det(l + A)) =-(trA)^ + 0(A='). (96) 
Next, note that for small e we have 

(1 + £)-V3 ^ I _ 1/3^ + 2/95' + 0{e^). (97) 
Combining this with (p^ . we get 

|(tr(TTA) - 3) = |((det(l + £:^))-^'hi{l + A) - 3 

= |((1 - l/3trA + l/18(trA)2 + l/6tr(A2) + 0(A=^))(3 + trA) - 3 

= ^(tr(A2) - l/3(trA)2) + 0{A^) = ^tr((A°)2) + 0(A=^). (9^ 

Finally, (^5]) follows from (1551) . using (1^ and 



Appendix B 

Let A, B G M and || A — B|| — 0. Let us prove, for instance, that 

B-^ : (A-B) = 0(||A-Bf). (99) 
Indeed, since det(-) is a smooth function, we have 

det(A) = det(B) + ; (A - B) + 0(11 A - Bf). (100) 

aB 

Next, using the Jacobi formula, we get 

det(A) = det(B) + det(B)B'T : (A - B) + 0(||A - Bf ). (101) 
Finally, taking into account that det(A) = det(B) = 1 and B""'" = B~^, we obtain (|99|) . 
Remark. Note that for the tangential space T-qM to the manifold M in Sym we have 

TbM = {X e Sym | B"^ : X = 0}. (102) 
Thus, relation (El) implies that lim ((A - B)/||A - B||) G TrM (if the limit exists). 

Appendix C 

We need to construct a function q{6) such that for small A 

q{e) > <l>(Cp\ Cp^) + 0(A) for all (Cp\ cf^) G 3(6, A). (103) 
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Let (Cp\ Cp^) G S{e, A). It follows from Appendix B, that 



(104) 



By X denote the orthogonal projection of C[^'' ^ — Cp"* ^ on the tangential space T^w-iM. 
Using f ll04p . we get for X 



X = C 



_ C(2)-1 



c 



_ . ^(1) 



cr 



ic 



{1)||2 



Moreover, since tr(xcS^^x) = ||X(CpV/2||2^ ^ave 

0(A3) 



(1) _ n(i)-i_n(2)-i + o(/N^2^)^ 

(105) 



trfXCp^X 



0(A). 



Substituting fllOSp in (168|) and taking fll06p into account, we obtain 



$(c«,cp)) 



-2(^4:^^:X)(l:X) 



tr 



(xcp^x 



+ 0(A). 



(106) 



(107) 



Thus, we define q{9) as 



q{e):= max g(C«), g(C(^)) 
||(c«-i)D||<e 



max 



tr 



(xcp)x 



(10^ 



The function q{C\^^) can be evaluated as follows. First, for each X introduce Y 



X(CPV/2_ Next, define a vector space T := {Y e Sym \ (CpV^^ : Y = 0}. Thus 



X G T^w-iM 



tr 



(xcp^x 



(i)-1ad 



where 



Bi := -2{C\ 



/p(l)-l\D 



Y G T, 



X = Bi : Y, 1 : X = B2 : Y, 



B2 := (C[ 



(l)\-l/2 



ll(cP"')°l, 

Therefore, 

[(b,:Y)(b,:Y 
Next, we compute the orthogonal projections of Bi and B2 on T: 

^{l)^l/2\fr^{l)^l/2_ ^ 



Bl :=B,- B,:(CP^)V^ (C| 



||(Cp))V^IP 



kG{l,2}. 



(109) 
(110) 

(111) 
(112) 

(113) 
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2x 10^ 
Ix 10^ 



0.01 0.02 0.03 

Fig. 6. Function q{0). 



Thus, 



max 

YeT,||Y||=l 



Y : sym(B? 



B :Y 



Amax(sym(B?(8)B° 



(114) 



where Ainax(sym(^B5 ® -^2)) maximal eigenvalue of the symmetric operator 

sym^B^^Bj) : Sym — > Sym. Obviously, the same maximal eigenvalue has its restriction 
on T° = Span{B°, B2}. It can be easily seen that 



sym(B; ® B") (B?) = 1/2(B? : B0)B? + 1/2(B? : BI)BI 



(115) 



sym(B? ® B°) (B^) = 1/2(B° : B0)B? + 1/2(B? : B^B^ 



(116) 



Therefore, the matrix of the restricted operator with respect to the basis {B^, B2} has 
the following form 



A :-- 



1 ( B; : BO B" : B" 
B;:B; B?:B^ 



;ii7) 



Both eigenvalues of A are real, since A represents a symmetric tensor. Finally, 



q{Cl'^) = A^ax(sym(B? ® B°)) = A^ax(^). 



(11^ 



Note that q{C\ ) is a continuous function of C- . Therefore, the maximum q{9) = 
max q(C^^^) is well defined. We compute it by the brutal force method. More- 

||((-;(1)-1)D||<5| 

over, the following parametrization can be used to simplify the computations. For any 
tensor Cp"* there exists a cartesian coordinate system and real numbers A\ A^ > such 
that the matrix of Cp'* takes the diagonal form diag(A^, A^, 1/(A^A^)). The function q{6) 
is plotted on the figure E] for <9 < 0.03. 
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